#ifndef CUBICBSPLINES_H
#define CUBICBSPLINES_H

#include "Spline.h"
#include "B_Basis_function.h"
//完全样条

class CubicBSpline_complete:public BSplines{
    public:
    std::vector<double>knots;
    std::vector<double>f;
    std::vector<double>coefficients;
    double m_1;
    double m_n;

    public:
    CubicBSpline_complete(std::vector<double>knots,std::vector<double>f,double m_1,double m_n);
    void CalculateCoefficients() ;
    double B_value(double x_value);
    };

    CubicBSpline_complete::CubicBSpline_complete(std::vector<double>knots,std::vector<double>f,double m_1,double m_n):BSplines(knots,f),knots(knots),f(f),m_1(m_1),m_n(m_n){        
        
        coefficients.resize(knots.size()+2); 
        CalculateCoefficients();
    }
void CubicBSpline_complete::CalculateCoefficients(){
    std::cout<<"complete"<<endl;
    int n=knots.size();
    Eigen::MatrixXd A=Eigen::MatrixXd::Zero(n+2,n+2);
    Eigen::VectorXd b=Eigen::VectorXd::Zero(n+2);
    //三次样条，每个区间上都受到三个函数的影响
    for(int i = 1; i <= n; i++){
        for(int j = i-1; j <= i+1; j++){
            B_Basis_function B(3, j-2, knots);
            A(i, j) = B.value(knots[i-1]);
        }
        b[i] = f[i-1];
    }
    
    //边界条件 
    B_Basis_function B__1(3,-2,knots);
    B_Basis_function B_0(3,-1,knots);
    B_Basis_function B_1(3,0,knots);

    B_Basis_function B_n__2(3,n-3,knots);
    B_Basis_function B_n__1(3,n-2,knots);
    B_Basis_function B_n(3,n-1,knots);



    A(0,0)=B__1.m(knots[0]);
    A(0,1)=B_0.m(knots[0]);
    A(0,2)=B_1.m(knots[0]);
    A(n+1,n-1)=B_n__2.m(knots[n-1]);
    A(n+1,n)=B_n__1.m(knots[n-1]);
    A(n+1,n+1)=B_n.m(knots[n-1]);

    b[0]=m_1;
    b[n+1]=m_n;

    Eigen::VectorXd a=Eigen::VectorXd::Zero(n+2);
    a=A.colPivHouseholderQr().solve(b);
    //std::cout<<"A："<<endl<<A<<endl;
    // 在求解完线性方程组后，将系数存入
    coefficients.resize(n+2);
    for (int i = 0; i <=n+1; i++) {
        coefficients[i]=a[i];
       // std::cout<<"系数："<<coefficients[i]<<endl;
    }
}

double CubicBSpline_complete::B_value(double x_value){
        double result=0.0;

        int n=knots.size();
        for(int i=0;i<n+2;i++){
            B_Basis_function B(3,i-2,knots);
          //  std::cout<<"B_value: "<<B.value(x_value)<<std::endl;
            result += coefficients[i]*B.value(x_value);
           // std::cout<<"coefficients: "<<coefficients[i]<<std::endl;
            
        }
       // std::cout<<"result: "<<result<<std::endl;
        return result;
}


    


//---------------------------自然样条------------------------------------------//


//自然样条
class CubicBSpline_natural:public BSplines{
    public:
    std::vector<double> knots;
    std::vector<double>f;
    std::vector<double>coefficients;

    public:
    CubicBSpline_natural(std::vector<double>knots,std::vector<double>f);
    void CalculateCoefficients();
    double B_value(double x_value);
    };

    CubicBSpline_natural::CubicBSpline_natural(std::vector<double>knots,std::vector<double>f):BSplines(knots,f),knots(knots),f(f){
        coefficients.resize(knots.size()+2); 
        CalculateCoefficients();
        
    }
void CubicBSpline_natural::CalculateCoefficients(){
    int n=knots.size();
    std::cout<<"natural"<<endl;
    Eigen::MatrixXd A=Eigen::MatrixXd::Zero(n+2,n+2);
    Eigen::VectorXd b=Eigen::VectorXd::Zero(n+2);
   for(int i = 1; i <= n; i++){
        for(int j = i-1; j <= i+1; j++){
            B_Basis_function B(3, j-2, knots);
            A(i, j) = B.value(knots[i-1]);
        }
        b[i] = f[i-1];
    }
    //边界条件 
    B_Basis_function B__1(3,-2,knots);
    B_Basis_function B_0(3,-1,knots);
    B_Basis_function B_1(3,0,knots);

    B_Basis_function B_n__2(3,n-3,knots);
    B_Basis_function B_n__1(3,n-2,knots);
    B_Basis_function B_n(3,n-1,knots);

    A(0,0)=B__1.M(knots[0]);
    A(0,1)=B_0.M(knots[0]);
    A(0,2)=B_1.M(knots[0]);
    A(n+1,n-1)=B_n__2.M(knots[n-1]);
    A(n+1,n)=B_n__1.M(knots[n-1]);
    A(n+1,n+1)=B_n.M(knots[n-1]);

    b[0]=0;
    b[n+1]=0;
  // std::cout<<"A："<<A<<endl;
    Eigen::VectorXd a=Eigen::VectorXd::Zero(n+2);
    a=A.colPivHouseholderQr().solve(b);
    // 在求解完线性方程组后，将系数存入
   // std::cout<<"系数："<<a<<endl;
    for (int i = 0; i <=n+1; i++) {
        coefficients[i]=a[i];
    }
    }
double CubicBSpline_natural::B_value(double x_value){
        double result=0.0;

        int n=knots.size();
        for(int i=0;i<n+2;i++){
            B_Basis_function B(3,i-2,knots);
          //  std::cout<<"B_value: "<<B.value(x_value)<<std::endl;
            result += coefficients[i]*B.value(x_value);
            //std::cout<<"coefficients: "<<coefficients[i]<<std::endl;
            
        }
        //std::cout<<"result: "<<result<<std::endl;
        return result;
}

//---------------------------周期样条------------------------------------------//


//周期样条
class CubicBSpline_periodic:public BSplines{
public:
    std::vector<double> knots;
    std::vector<double> f;
    std::vector<double> coefficients;

public:
    CubicBSpline_periodic(std::vector<double> knots,std::vector<double> f);
    void CalculateCoefficients();
    double B_value(double x_value);
};
CubicBSpline_periodic::CubicBSpline_periodic(std::vector<double> knots,std::vector<double> f):BSplines(knots,f),knots(knots),f(f){
    coefficients.resize(knots.size()+2); 
    CalculateCoefficients();
    }

void CubicBSpline_periodic::CalculateCoefficients() {
    std::cout << "periodic" << std::endl;

    int n = knots.size();

    // Check if f(a) == f(b)
    if (std::abs(f[0] - f[n-1]) > 1e-10) {
        std::cout << "f(a) must = f(b)!" << std::endl;
        return;
    }

    // Construct the periodic system of equations
    Eigen::MatrixXd A = Eigen::MatrixXd::Zero(n + 4, n + 2);
    Eigen::VectorXd b = Eigen::VectorXd::Zero(n + 4);

    // Internal points equations
    for (int i = 2; i <= n + 1; i++) {
        for (int j = i - 2; j <= i; j++) {
                B_Basis_function B(3, j-2, knots);
                A(i, j) = B.value(knots[i - 2]);
                //std::cout << "A[" << i << "][" << j << "] = " << A(i, j) << std::endl;
            }
        b(i) = f[i - 2];
    }

    // Periodicity constraints
    B_Basis_function B__1(3,-2,knots);
    B_Basis_function B_0(3,-1,knots);
    B_Basis_function B_1(3,0,knots);

    B_Basis_function B_n__2(3,n-3,knots);
    B_Basis_function B_n__1(3,n-2,knots);
    B_Basis_function B_n(3,n-1,knots);

    // First derivative periodicity
    A(0, 0) = B__1.m(knots[0]);
    A(0, 1) = B_0.m(knots[0]);
    A(0, 2) = B_1.m(knots[0]);

    A(n + 2, n - 1) = B_n__2.m(knots[n - 1]);
    A(n + 2, n ) = B_n__1.m(knots[n - 1]);
    A(n + 2, n+1) = B_n.m(knots[n - 1]);

    b(0) = b(n + 2);

    // Second derivative periodicity
    A(1, 0) = B__1.M(knots[0]);
    A(1, 1) = B_0.M(knots[0]);
    A(1, 2) = B_1.M(knots[0]);

    A(n + 3, n - 1) = B_n__2.M(knots[n - 1]);
    A(n + 3, n ) = B_n__1.M(knots[n - 1]);
    A(n + 3, n+1) = B_n.M(knots[n - 1]);

    b(1) = b(n + 3);

  //std::cout << "A:\n" << A << std::endl;
   // std::cout << "b:\n" << b << std::endl;

    // Solve the linear system
    Eigen::VectorXd x = A.colPivHouseholderQr().solve(b);
   //std::cout << "Coefficients:\n" << x << std::endl;

    // Store coefficients safely
    coefficients.resize(n + 2); // Ensure sufficient size
    for (int i = 0; i < n + 2; i++) {
        coefficients[i] = x[i];
    }
}
double CubicBSpline_periodic::B_value(double x_value){
        double result=0.0;

        int n=knots.size();
        for(int i=0;i<n+2;i++){
            B_Basis_function B(3,i-2,knots);
           // std::cout<<"B_value: "<<B.value(x_value)<<std::endl;
            result += coefficients[i]*B.value(x_value);
            //std::cout<<"coefficients: "<<coefficients[i]<<std::endl;
            
        }
      //  std::cout<<"result: "<<result<<std::endl;
        return result;
}
#endif